Matlab Codes

Chapter 3

Matlab code 3.2: Matlab file “Figure 3-9.m”

%--------------------------------------------------------------------

% This code can be used to generate Figure 3.9-3.11

%--------------------------------------------------------------------

clear all;

close all;

%% the theoretical curve of the vibrating target micro-Doppler characteristic(LFM radar)

c = 3e8;

j = sqrt(-1);

fc = 10e9; % carrier frequency of transmitted signal

cord = 1000*[3 4 5]; % coordinates of local coordinate system's origin in the radar coordinate system

colo = [0 0 0]; % coordinates of radar in the radar coordinate system

R0 = cord-colo;

alpha = 0.9273; % azimuth angle

% alpha = atan(cora(2)/cora(1));

beita = 0.7854; % elevation angle

% beita = atan(cora(3)/(sqrt(cora(1)^2+cora(2)^2)));

alpha_p = pi/5; % azimuth angle of vibrating axis

beita_p = pi/4; % elevation angle of vibrating axis

v = 0; % translational velocity of target

dv = 0.1; % vibration amplitude

fv = 3; % vibration frequency

prf = 1000; % pulse repetition frequency

pri = 1/prf; % pulse repetition interval

t = 1; % radar illumimated time

dt = 0:pri:t-pri; % time sampling interval

temp = cos(alpha-alpha_p)*cos(beita)*cos(beita_p)+sin(beita)*sin(beita_p);

fmd = (4*pi*fc*fv*dv/c)*temp*cos(2*pi*fv*dt); % theoretical micro-Doppler frequency

%% the simulation result of the vibrating target micro-Doppler characteristic(LFM radar)

r = zeros(length(colo),length(dt)); % distance between the scatterers and radar

m = length(dt);

for i = 1:m

r(:,i) = R0'+dv*sin(2*pi*fv*dt(i))*[cos(alpha_p)*cos(beita_p);sin(alpha_p)*cos(beita_p);sin(beita_p)];

end

r_abs = sqrt(sum(r.^2,1));

s = exp(j*2*pi*fc*2*r_abs'/c); % echo signal matrix

N = 200; % number of Gabor coefficients in time

Q = 100; % degree of oversampling

tfr1 = tfrgabor(s.',N,Q); % Gabor transform

tfr_r1 = fftshift(tfr1,1);

%% range-slow time profile

B = 500e6; % bandwidth

tp = 5e-6; % pulse duration

mu = B/tp; % frequency modulation rate

xmin = sqrt(sum((cord-colo).^2))-10*dv; % the minmum detectable distance

xmax = sqrt(sum((cord-colo).^2))+10*dv; % the maximum detectable distance

ts = 2*xmin/c; % the starting point of the time sampling

te = 2*xmax/c+tp; % the ending point of the time sampling

fs = B; % fast time sampling frequency

kdt = 1/fs; % fast time sampling interval

nt = 2*ceil((te-ts)/(2*kdt)); % number of fast time sampling

df = fs/nt; % frequency sampling interval

tk = ts+(0:nt-1)*kdt; % fast time

rw = zeros(length(colo),length(dt));

for i = 1:m

rw(:,i) = sqrt((R0'+dv*sin(2*pi*fv*dt(i))*[cos(alpha_p)*cos(beita_p);sin(alpha_p)*cos(beita_p);sin(beita_p)])'...

*(R0'+dv*sin(2*pi*fv*dt(i))*[cos(alpha_p)*cos(beita_p);sin(alpha_p)*cos(beita_p);sin(beita_p)]));

end

for i = 1:length(colo)

td = tk(:)*ones(1,m)-ones(length(tk),1)*rw(i,:)*2/c; % time of arrival

end

s = exp(j*2*pi*fc*td+j*pi*mu*td.^2);

R0_abs = sqrt(sum(R0.^2));

td0 = tk(:)*ones(1,m)-2*R0_abs/c;

s0 = exp(j*2*pi*fc*td0+j*pi*mu*td0.^2); % reference signal

sdt = s.*conj(s0); % dechirp

sdf = fftshift(fft(sdt),1);

figure (1)

x = dt;

y = -(df*(0:nt-1)-B/2)*c/(2*mu); % turn the fast time frequency-slow time plane to the range-slow-time plane

contour(x,y,abs(sdf))

xlabel('Slow time (s)')

ylabel('Range (m)')

axis([0,1,-8,8])

%% extraction of frequency domain

N = 200;

Q = 200;

tfr2 = tfrgabor(sdf(nt/2+1,:)',N,Q);

tfr_r2 = fftshift(tfr2,1);

figure (2)

contour(linspace(0,t,length(tfr_r2(1,:))),linspace(-50,50,length(tfr_r2(:,1))),tfr_r2,30)

xlabel('Slow time (s)')

ylabel('Frequency (Hz)')

axis([0,1,-40,40])

%% extraction of time domain

N = 200;

Q = 100;

tfr3 = tfrgabor(sdt(nt/2+1,:)',N,Q); % Gabor transform

tfr_r3 = fftshift(tfr3,1);

figure (3)

contour(linspace(0,t,length(tfr_r3(1,:))),linspace(-50,50,length(tfr_r3(:,1))),tfr_r3,30)

xlabel('Slow time (s)')

ylabel('Frequency (Hz)')

axis([0,1,-40,40])